subroutine r2abij_t2aeij_fbe(r2, t2, fae, & 
               icore, fact)
!
use mod_size
use mod_ioff
use mod_orbit
use mod_iop 
implicit none
!
real*8, intent(in) :: fact
real*4, intent(in) :: r2(*) 
real*4, intent(in) :: fae(*)
!
real*4, intent(inout) :: t2(*), icore(*)
!
integer :: i, j, a, b, nsize
integer :: i0, i1, i2, i3
integer, external :: idsymsoc
real*8, external :: dnrm2, nrm2so
!this is to do r2(ab,ij) = r2(abij) + P(ab)t2(ae,ij)*f(be)
i0 = 1               
i1 = i0 + nvvoo*isd 
i2 = i1 + nvvoo*isd
i3 = i2 + nvvoo*isd
!
!AAAA===
!t2(AE,IJ)*Fae(BE)
!
call itranspso(t2(ioi3(1)), icore(i1), vrta, vrta, 0, popa, popa, 0, 1) 
call iexpso(icore(i1), icore(i0), popa, popa, 0, 0, vrta, vrta, 0, 1, 1) 
call VpqrmVsm_to_Vpqrs(icore(i0), fae(ifvv(1)), 1, 1, vrta, & 
             popa, popa, 0, vrta, vrta, icore(i1), 1, 1.d0, 0.d0, 0)
!--tmp(i<j,a<b) = icore(i<j,ab)-icore(i<j,ba)
nsize = idsymsoc(1, vrta, vrta, 1, popa, popa, 0)
!write(6,*)'norm of icore(i1)', nrm2so(nsize, icore(i1), 1)
call iassymso(icore(i1), icore(i0), 1, popa, popa, vrta, vrta, 0)
nsize = idsymsoc(1, vrta, vrta, 0, popa, popa, 0)
!write(6,*)'norm of icore(i2)', nrm2so(nsize, icore(i0), 1)
call itranspso(icore(i0), icore(i1), popa, popa, 0, vrta, vrta, 0, 1) 
!write(6,*)'norm of icore(i3)', nrm2so(nsize, icore(i1), 1)
!write(6,*)'norm of r2 before', fact, nrm2so(nsize, r2(ioi3(1)), 1)
call axpyso(nioi31, fact, icore(i1), 1, r2(ioi3(1)), 1)
!write(6,*)'norm of r2 after', nrm2so(nsize, r2(ioi3(1)), 1)

!
!ABAB===
! t2(Ae,Ij)*Fae(be) - t2(bE,Ij)*Fae(AE)
!=t2(Ae,Ij)*Fae(be) + t2(Eb,Ij)*Fae(AE)
call itranspso(t2(ioi3(2)), icore(i0), vrta, vrtb, 1, popa, popb, 1, 1) 
call VpqrmVsm_to_Vpqrs(icore(i0), fae(ifvv(2)), 1, 1, vrtb, & 
                 popa, popb, 1, vrta, vrtb, icore(i1), 1, 1.d0, 0.d0, 1) 
!
call isymtrso(icore(i0), icore(i2), popa, popb, 1, 0, vrta, vrtb, 1, 1, 1)
call VpqrmVsm_to_Vpqrs(icore(i2), fae(ifvv(1)), 1, 1, vrta, & 
                 popa, popb, 1, vrtb, vrta, icore(i0), 1, 1.d0, 0.d0, 1) 
nsize = idsymsoc(1, vrta, vrtb, 1, popa, popb, 1) 
!write(6,*)'norm of icore(i1)', nrm2so(nsize, icore(i0), 1)
call isymtrso(icore(i0), icore(i2), popa, popb, 1, 0, vrtb, vrta, 1, 1, 1)
call axpyso(nioi32, 1.d0, icore(i2), 1, icore(i1), 1) 
!
call itranspso(icore(i1), icore(i0), popa, popb, 1, vrta, vrtb, 1, 1) 
call axpyso(nioi32, fact, icore(i0), 1, r2(ioi3(2)), 1) 
!
!BBBB===
!t2(ae,ij)*Fae(be)
!
call itranspso(t2(ioi3(3)), icore(i1), vrtb, vrtb, 0, popb, popb, 0, 1)
call iexpso(icore(i1), icore(i0), popb, popb, 0, 0, vrtb, vrtb, 0, 1, 1)
call VpqrmVsm_to_Vpqrs(icore(i0), fae(ifvv(2)), 1, 1, vrtb, &
             popb, popb, 0, vrtb, vrtb, icore(i1), 1, 1.d0, 0.d0, 0)
nsize = idsymsoc(1, vrtb, vrtb, 1, popb, popb, 0) 
!write(6,*)'norm of icore(i1)', nrm2so(nsize, icore(i1), 1)
!--tmp(i<j,a<b) = icore(i<j,ab)-icore(i<j,ba)
call iassymso(icore(i1), icore(i0), 1, popb, popb, vrtb, vrtb, 0)
call itranspso(icore(i0), icore(i1), popb, popb, 0, vrtb, vrtb, 0, 1)
call axpyso(nioi33, fact, icore(i1), 1, r2(ioi3(3)), 1)
!
return
end
